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The  Independent  block  structure  of  regenerative  processes  and  the  knoUU  l!uu“ 
vergence  rates  of  the  means  of  ratio  estimators  are  exploited  to  produce  bias-free 
regression-adjusted  estimates  (AOAe's)  for  regenerative  simulations.  Direct  assessments 
of  the  variances  of  the  estimates  are  obtained,  as  well  as  indications — both  formal  and 
graphical — of  normality  and  symmetry,  or  lack  of  it,  in  the  distribution  of  the  estimates. 


1.  INTRODUCTION 

In  regenerative  simulations  one  exploits  the  fact  that  the  sample  path  of  the 

simulated  process,  say  {W^  i  >  0},  can  be  divided  into  independent  and  identically 

distributed  (i.i.d.)  blocks  of  lengths  {x^,  j  _>  1} .  For  example,  let  the  quantity  of 

interest  be  the  expected  stationary  waiting  time,  E(W),  in  a  queue  and  let  Y^  be  the 

sum  of  the  waiting  times  observed  in  the  jth  block.  Then  the  usual  regenerative 

estimate,  Ae(n) ,  of  E(W),  formed  from  n  blocks  is  Ae(n)  •  Y/x  where  Y  and  x 

are  the  averages  of  {Y. .  Y  } ,  and  {x,,...,  x  },  respectively.  This  estimate  con- 

l  n  in 

verges  to  E(W)  as  n  gets  large  and  the  so-called  "statistics"  of  regenerative 
simulation  exploit  the  fact  that  4e(n)  is  asymptotically  unbiased  and  normally  dis¬ 
tributed  with  a  variance  of  a  known  form  which  can  be  estimated  from  the  data  [Crane, 
and  Iglehart  (1975);  Crane  and  Lemoine  (1977);  Iglehart  (1978)]. 

An  advantage  of  the  regenerative  estimate  is  that  it  eliminates  the  problem  which 
arises  in  straightforward  simulations  of  choosing  initial  conditions  so  as  to  make 
(W  ,  i  >  0}  stationary  and  W,  the  sample  path  average,  unbiased.  It  is  crucial  to 
understand  that  the  initial  conditions  are  of  importance  only  for  finite  sample 

*Work  supported  in  part  by  the  Office  of  Naval  Research  under  Grant  NR-42-84. 
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sises  so  that  the  advantage  of  Ae(n)  over  W  disappears  asymptotically.  Given 
this  advantage,  the  regenerative  estimate  is  usually  employed  and  a  single  sample  path 
is  used  for  as  long  as  one  needs  to  simulate  in  order  to  achieve  a  given  precision.  At 
that  point  the  simulator  assumes  that  the  point  estimate  Ae(n)  is  both  unbiased  and 
normally  distributed.  There  are  several  problems  involved  in  this  "flying-blind" 
procedure  besides  the  fact  that  the  assumptions  of  normality  and  no  bias  are  typically 
not  verified  and  may  not  be  true.  Most  notably,  it  is  not  known  how  to  adjust  confi¬ 
dence  interval  widths  to  allow  for  the  fact  that  the  variance  is  estimated  from  the  data 
and  is  usually  highly  positively  correlated  with  the  point  estimate.  This  correlation 
is  particularly  troublesome  in  sequential  procedures  which  use  a  relative  confidence- 
interval-width  criterion  for  determining  simulation  run  lengths. 

A  method  which  is  preferable  to  using  a  single  sample  path  would  be  to  generate 
m  l.i.d.  sample  paths,  form  an  averaged  estimate  from  the  m  estimates  obtained  and 
then  produce  the  usual  standard  deviation  estimate  for  a  sample  mean.  Unfortunately 
taking  m  sample  paths  of  length  n/m  only  aggravates  the  concern  with  bias  and 
normality  of  the  estimate  which  occurs  with  both  the  straightforward  and  regenerative 
methods. 

In  this  paper  we  exploit  two  aspects  of  the  regenerative  structure  of  the  simu¬ 
lated  process  to  resolve  this  dilemma.  The  first  is  that  the  regenerative  blocks 
are  l.i.d.  so  that  a  single  simulation  of  n  cycles  can  always  be  sectioned  into  k 
i.l.d.  simulations  of  n^  ■  n/k  blocks  per  section  for  multiple  values  of  k.  Thus 
from  a  single  simulation  of  n  blocks,  unbiased  (but  correlated)  estimates  of 
E(A£(n^))  can  be  obtained  for  the  multiple  values  of  n^.  Secondly,  unlike  the 
estimate  W,  the  regenerative  estimate  has  known  bias  structure,  l.e., 

E(Ae(n))  ■  E(W)  +  g^/n  +  g0/n^+  •••  .  Combining  these  two  facts  and  using  regression 
techniques  an  estimate  and  a  graphical  picture  of  the  bias  structure  can  be  obtained. 
More  Importantly,  a  regression  adjusted  regenerative  estimate.  AdAe(n),  for  E(W)) 

can  be  produced  which  corrects  for  the  bias. 
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This  procedure  can  itself  be  replicated,  say,  m  times .producing  m  i.i.d. 

estimates  A4Ae(l,n) . m,n)  with  low  bias  and  the  value  E(W)  can  then  be 

estimated  by  (VuV ie(m,n),  the  average  of  these  m  ,n)  's.  Furthermore  plots  showing 

the  evolution  to  the  unbiased  state  can  be  produced  along  with  box  plots  displaying  the 
distributions  of  the  estimates.  Running  normal  plots  with  normal  test  statistics  can 

*  also  be  given  for  various  sub run  lengths  n^.  These  plots  allow  the  user  to  judge  the 
sub8ample  length  n  at  which  Aa/iz(l,n)  , . . .  ,AoAe(m,n)  can  be  assumed  to  be  normally 

w 

distributed.  Using  this  sample  of  m  Aa/iz's  a  standard  deviation  estimate,  6tiOAt( m,n), 
for  the  averaged  KCUite,  £UutAZ(m,n) ,  is  obtained  from  which  a  confidence  Interval  for 

the  steady  state  parameter  is  derived.  Furthermore  the  sample  of  m,  /UZA.e(l,n) . 

A/Vlz( m,n),  is  available  to  do  comparisons  and  rankings  and  to  form  multidimensional  con¬ 
fidence  regions  (when  more  than  one  parameter  is  of  interest)  using  standard  normal 
theory  statistical  methodology  or  nonparametric  methodology. 

Other  techniques  which  reduce  bias  in  regenerative  estimates  are  the  jackknife 

and  Tin  estimates  [Miller  (1970);  Iglehart  (1975)]  and  estimates  based  on  renewal 

theoretic  properties  [Heldelberger  (1978),  and  Meketon  (1979)].  As  generally  applied 

these  techniques  remove  only  the  1/n  term  in  the  bias  expansion  (the  jackknife  can 
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be  modified  to  remove  the  1/n  term  at  a  great  computational  cost  and  an  uncertain 
effect  on  the  variance  of  the  resulting  estimate).  However,  no  guidance  is  given  as 
to  when  this  bias  reduction  is  sufficient  and  the  jackknife  is  known  to  inflate  the 
variance  in  finite  samples  [Efron  and  Stein  (1978);  Goodman,  Lewis  and  Robbins  (1971)]. 
Furthermore,  the  problems  of  correlation  between  the  point  and  variance  estimates  and 
of  detecting  normality  are  not  relieved  by  these  alternative  estimators.  The  regres¬ 
sion  and  graphical  procedures  given  in  Sections  3  and  4  could  be  applied  to  estimates 

*  for  which  the  1/n  term  has  already  been  removed,  by,  for  instance,  jackknifing,  al¬ 
though  this  seems  redundant.  In  practice  the  bias  becomes  insigif leant  long  before  the 
assumption  of  normality  is  valid,  particularly  when  bias  reducing  techniques  are  applied. 

The  work  in  this  paper  is  related  to  that  of  Fishman  (1977)  who  used  replica¬ 
tions  with  the  Tin  estimate  to  reduce  bias,  and  a  normality  test.  The  strength  of  the 


3 


regression  adjusted  regenerative  estimate  in  conjunction  with  the  graphics  is  that  one 
obtains  s  very  clear  picture  of  how  the  bias  is  changing  and  at  what  point  it  is  insig¬ 
nificant.  In  addition,  the  running  normal  plots  show  the  convergence  to  normality. 

To  apply  the  methodology  given  in  this  paper  in  practice  requires  a  relatively 
sophisticated  user  to  interactively  interpret  the  output.  This  sharply  contrasts  with 
previous  automated  sequential  procedures  [Fishman  (1977)  and  Lavenberg  and  Sauer  (1977)], 
which  require  no  human  intervention,  but  which  provide  the  user  with  far  less  informa¬ 
tion,  particularly  about  the  dynamics  of  the  simulation. 

2.  OBJECTIVES  AND  OUTLINE 

The  various  regression  adjusted  estimates  are  described  in  Section  3  of  the  paper 
and  the  graphics  which  has  been  developed  to  go  with  them  are  described  in  Section  4. 

In  Section  5  we  describe  a  protocol  which  directs  the  simulator  in  the  use  of  the  above 
methodology. 

The  objective  of  this  protocol  is  to  produce,  for  a  given  precision  or  total 
available  computing,  estimates  kake( j,m),  j  “  l,...,m,  from  the  shortest  possible 
number  of  blocks  per  replication  such  that  the  AaA&(j ,m) 's  are  approximately 
(1)  unbiased; 

(11)  normally  distributed. 

The  minimisation  of  the  number  of  blocks  in  each  replication  is  done  in  order  to  optimize 
the  number,  m,  of  replications  of  the  simulation  from  which 

(i)  a  standard  deviation  estimate  4kake(m,n)  for  the  averaged  stake  estimates, 
Okake(m,n) ,  is  obtained; 

(11)  a  confidence  interval  for  the  parameter  is  obtained  using  4kOke(n,n) 
and  0k0ke(m,n) ; 

(ill)  a  sample  of  m  MU te(j,n)'s  is  obtained  to  do  comparisons,  rankings,  etc.  using 
standard  normal  theory  statistical  methodology  or  nonparametric  techniques. 

Of  course  in  some  situations  it  may  well  be  that  not  enough  independent  blocks  can 
be  simulated  to  be  able  to  obtain  replications.  The  utility  of  the  proposed  procedure 
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then  Is  that  even  though  the  user  is  forced  back  onto  the  regular  regenerative  estimate 
A£(n)  he  is  not  doing  so  blindly.  Moreover,  the  single  regression  adjusted  regenerative 
estimate  >uv t£(n)  provides  him  with  an  essentially  bias  free  alternative  to  the  regen¬ 
erative  estimate. 

In  Section  6  examples  of  the  use  of  the  protocol  and  its  graphical  displays  are 
given  for  three  queueing  models  and  in  Section  7  simulation  results  on  the  relative 
properties  of  A.£(nm)  and  a/iasiz(m,n)  are  presented.  Section  8  contains  a  summary. 
Indicates  related  topics  for  further  research,  and  gives  examples  of  situations  other 
than  regenerative  simulations  for  which  this  methodology  applies.  Table  4  gives  a 

summary  of  the  notation  used  in  the  paper. 

3.  THE  REGRESSION  ADJUSTED  REGENERATIVE  ESTIMATE 

Under  broad  conditions  the  expected  value  E(A£(n))  has  the  form 

E(A£(n))  -  80  +  Sjn-1  +  ^n"2  +•••,  (3.1) 


where  8^  *»  E(w) .  Let  a  section  of  n  regenerative  blocks  be  broken  up  into  subsections 
of  lengths  n^  for  k  -  1, ...,K  and  let  m^  «  [n/n^],  where  [x]  denotes  the  greatest 
integer  less  than  or  equal  to  x.  For  each  k,  m^  regenerative  estimates,  -te(j,n^), 

J  ■  l,...,!!^,  are  formed  from  disjoint  blocks  of  n^  cycles. 

Regression  adjusted  regenerative  estimates  are  formed  by  estimating  Bq  in 

(3.1)  using  the  AeO.n^’s  for  j  -  1 . n^,  k  -  1,...,  K,  and  the  equation  (3.1) 

to  some  degree  d  ■  1,2,  or  3.  This  can  be  done  in  three  distinct  ways: 

(i)  By  forming  average  regenerative  estimates,  aA-e^.n^) ,  for  k  »  1,...,K 

m. 

where  OieCn^,!^)  -  *£(3,1^)/:^,  and  using  these  as  dependent  variables  in  the 

regression  on  6^  +  +•••+  Bdn^d.  The  averaging  will  tend  to  make  the  dependent 

variables  more  normal  than  the  individual  A£(j,nji)'s  and  since  var(A£(j ,n^)) 

+  d/n^  +•••»  then  we  neglect  the  higher-order  terms  var(aA.Z(u^tt^)  is 

approximately  y/n  for  all  k  *  1 . K.  Thus  an  unweighted  least  squares  estimate 

of  6q  will  produce  $q,  the  40*£(n)  estimate  unbiased  out  to  terms  up  to  1/n^. 


(ii)  By  using  the  >te(j  .n^/m^  for  j  ■  1 . m^,  k  ■  1, . . .  ,K  in  a  weighted 

regression  on  Bq/oi^  +  (B^/m^n^*  +•••+  (B^/m^n^.  ^lis  *8  e<luivalent  to  (D  if 
least  squares  estimation  of  Bq,...,  3^  is  performed. 

(Hi)  Since  it  will  be  seen  in  later  sections  that  4&(j,n^)  is  highly  non-normal 
for  small  n^,  robust  regression  methods  could  be  used  to  eliminate  the  influence  of 
outliers.  However,  since  the  4£(j,n^)'s  are  generally  not  symmetric  random  variables, 

A 

the  8q  estimate  will  generally  be  biased. 

In  this  paper  the  KOHi.  estimate  is  always  Bq  obtained  as  in  (i) ,  which  like 

-te(n)  is  a  statistic  involving  (Y, ,...,Y  )  and  (x, . t  ).  The  use  of  robust 

in  in 

regression  techniques  will  be  described  elsewhere. 

There  are  still  three  problems  which  make  the  regression  non-standard: 

(i)  How  large  should  K  be,  given  that  it  must  be  at  least  greater  than  d+1? 

(ii)  The  o4.e(c^,n^)  are  correlated,  which  seems  to  mitigate  against  making  K 
too  large. 

(iii)  Given  that  K  is  chosen,  how  should  the  n^  be  chosen?  This  is  an  experimental 
design  problem. 


These  are  difficult  and  comprehensive  analytic  problems,  not  considered  here. 
However,  it  is  possible  to  show  that  the  correlation  between  <ZA.6(m^,n^)  and 
a/tfc(mj  .n^)  converges  to  one  if  n  ♦  00  in  such  a  way  that  m^  and  m^  are  fixed 
divisors  of  n.  Furthermore  in  this  case  var('UM.e(n))/varOte(n))  also  converges  to 
one.  For  small  to  medium  sample  sizes  one  suspects  that  the  correlation  will  inflate 
the  variance  of  the  KOK<L( n)  estimate  as  compared  to  the  variance  of  the  /te(n) 
estimate.  This  frequently  happens  in  small  samples  with  the  jackknife.  However , 
simulation  studies  presented  in  Section  7  show  that  var(/uW.e(n))  is  only  very 
slightly  larger  than  var(*e(n))  and  that  this  result  is  robust  with  respect  to  the 
choice  of  K  and  the  i^'s.  Of  course,  choosing  too  small  will  result  in  non¬ 
normality  and,  possibly,  bias  (the  degree,  d,  being  inadequate). 

We  emphasize  here  that  the  preference  for  the  *a*e(n)  estimate  over  and 
above  a  jackknife  estimate  or  a  Tin  estimate,  is  that  its  components  allow  one  to 


study  the  evolution  of  the  regenerative  estimate.  The  graphics  to  do  this  are  dis¬ 
cussed  in  the  next  section. 

Of  course  one  utility  of  the  jackknife  is  that  it  can  be  used  with  small  samples 
that  a  jackknifed  variance  estimate  is  available.  However,  this  last  advantage  is  not  an 
issue  here  since  a  standard  variance  estimate  is  available  for  regenerative  estimate 
A£(n);  furthermore,  once  a  sectioning  into  m  sections  is  achieved,  even  this  can  be 
dispensed  with.  A  regression-based  estimate  of  the  variance  of  Ae(n)  can  be  constructed 
but  is  not  discussed  here. 

4.  DESCRIPTION  OF  THE  GRAPHICS 

Three  different  types  of  graphs  are  used  to  display  and  analyze  the  output. 

The  first  two,  the  Basic  Graph  and  the  Retrenched  Graph,  are  quite  similar  and  are 
concerned  with  identifying  the  evolution  of  bias,  skewness  and  departures  from  normality 
(in  the  form  of  individual  outliers  or  outlying  sample  paths)  in  the  regenerative  and 
fiaAZ  estimates.  The  third  graph,  the  running  normal  graph,  specifically  displays  the 
convergence  to  normality  of  the  estimates.  These  three  types  of  graphs,  described 
in  this  section, all  deal  with  the  data  of  primary  concern,  namely,  the  /te,  (UiQ. ,  MAS., 
and  a/UVie.  estimates.  They  are  novel  and  are  adaptations  of  ideas  which  are  pertinent 
in  other  estimation  contexts. 

(Other  plots  that  may  be  of  interest  in  a  regenerative  simulation  include 
histograms  and  box  plots  of  {Y^,  k  D  and  k  1) .  Scatter  plots  of 

{ (Yfc* Tk)»  k  2.  ^  show  the  high  correlation  that  typically  exists  between  Yfc  and  t^. 
These  last  three  plots  are  useful  in  explaining  the  skewed  distribution  often  exhibited 
by  regenerative  estimates.  However,  since  we  feel  their  importance  is  secondary  and 
they  are  well  known  they  will  not  be  discussed  here.) 

We  proceed  by  first  describing  the  box  plot  which  is  a  primary  component  of 
the  Basic  and  Retrenched  Graphs.  A  description  and  interpretation  of  the  Basic  and 

Retrenched  Graphs  comes  next,  followed  by  a  description  of  normal  plots  and  the 
third  type  of  graph,  a  sequence  of  normal  plots  called  a  running  normal  graph. 

A  protocol  for  the  interactive  use  of  these  graphs  is  outlined  in  Section  5  and 
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examples  of  Its  use  Is  given  in  Section  6.  The  actual  plots  in  this  paper  were  pro- 
duced  on  a  Tektronix  4013  terminal  using  an  APL  graphics  package.  The  user  may  wish 
to  refer  to  Table  4  for  a  summary  of  the  notation  used  in  the  paper. 

4.1.  The  Box  Plot 

A  convenient  and  compact  technique  for  displaying  the  distribution  of  a  batch 
of  data  is  the  box  plot.  These  plots  were  Introduced  by  Tukey  (1977),  although  we 
prefer  the  form  adapted  by  McNeil  (1977).  The  box  plot  is  an  excellent  tool  for 
Identifying  skewness  and  asymmetry  in  the  data  as  well  as  for  detecting  outliers. 

The  plot  is  therefore  well  suited  to  analyzing  the  often  highly  skewed  and  sometimes 
wild  output  of  queueing  simulations  (see,  e.g.,  the  left-most  box  plot  of  Figure  4C 
at  n  ■  50;  the  160  regenerative  estimates  which  make  up  the  box  plot  are  highly 
positively  skewed  and  clearly  non-normal) . 

Figure  1  shows  a  sample  box  plot.  The  body  of  the  box  plot  displays  the  median  j 

and  the  lower-  and  upper-quar tiles  of  the  data.  Let  q  represent  the  (estimated)  ) 

J 

Interquartile  distance;  q  *  upper  quartile  minus  lower  quartile.  Those  data  points  between^ 

j 

the  lower  quartile  minus  1.5q  and  the  lower  quartile  minus  q  are  marked  by  light  circles, 

as  are  those  data  points  between  the  upper  quartile  plus  q  and  the  upper  quartile 
plus  1.5q.  Those  values  below  the  lower  quartile  minus  1.5q  or  above  the  upper  quartile 

plus  1.5q  are  marked  by  dark  circles  and  are  meant  to  be  indicated  as  outliers.  For 
normally  distributed  data,  approximately  5%  of  the  points  should  be  marked  by  light 
circles,  but  only  about  1  in  200  should  be  marked  by  dark  circles.  One  problem  with 
box  plots  is  that  they  suppress  sample  size  and  this  can  lead  to  erroneous  inferences; 
see  McGill,  Tukey  and  Larsen  (1978)  for  extended  box  plots  which  take  account  of  this. 

For  samples  of  size  less  than  9  the  box  plots  are  not  very  meaningful  and  the  data  points  • 
are  laid  out  vertically  by  magnitude  for  visual  inspection. 

m 

4.2.  The  Basic  and  Retrenched  Graphs 

The  two  primary  graphs  used  in  the  output  analysis  are  the  Basic  Graphic  and 
the  Retrenched  Graph.  The  Basic  Graph  provides,  in  the  first  stage  of  the  protocol. 


8 


sample:  box  plot 

PLOT  OF  50  POINTS 

•  __  UPPER  QUARTILE  +  1 . 5* INTERQUARTILE  DISTANCE 
o 

v  __  LARGEST  VALUE  s  UPPER  QUARTILE  +  INTERQUARTILE  DISTANCE 

j-L  __  UPPER  QUARTILE 

_  __  MEDIAN 

—J  __  LOWER  QUARTILE 

__  SMALLEST  VALUE  *  LOWER  QUARTILE  -  INTERQUARTILE  DISTANCE 
°  __  LOWER  QUARTILE  -  1 .5* INTERQUARTILE  DISTANCE 

t 

INTERQUARTILE  DISTANCE  =  UPPER  QUARTILE  -  LOWER  QUARTILE 


FIGURE  1.  Example  box  plot  for  fifty  data  points  from  a  regenerative 

simulation.  The  interquartile  distance  equals  the  estimated 
upper  quartile  minus  the  estimated  lower  quartile.  The  light 
circles  are  data  points  which  fall  between  the  largest  value 
less  than  or  equal  to  the  upper  quartile  plus  the  inter¬ 
quartile  distance  and  the  upper  quartile  plus  1.5  times  the 
interquartile  distance.  The  dark  circles  are  data  with  values 
above  this  latter  point.  Similarly  for  the  lower  part  of 
the  box  plot. 
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preliminary  estimates  of  the  section  length,  n,  and  the  degree  of  the  polynomial,  d, 
needed  in  the  regression  to  reduce  the  bias  to  an  inconsequential  level.  As  the  total 
run  length  increases,  the  Retrenched  Graph  is  used  to  adjust  these  quantities  and  to 
monitor  the  distributional  convergence  of  the  ^e's,  asit' s  and  -’UZAe's. 

The  Retrenched  Graph  displays  the  results  of  performing  m  i.i.d.  replications 
of  the  regression  as  described  in  Section  3.  The  Basic  Graph  is  essentially  a  special 
case  of  the  Retrenched  Graph  with  m  ■  1.  Examples  of  Basic  Graphs  are  Figures  3A, 

3B,  4A  and  5A  while  Figures  3C,  3D,  3F,  3G,  31,  4B,  4C,  4D,  4E,  4G  and  5B  are  all 

Retrenched  Graphs.  The  following  discussion  describes  the  features  of  the  Retrenched 
Graph;  it  applies  as  well  to  the  Basic  Graph  with  m  set  equal  to  one.  Following  the 
description  is  a  discussion  on  interpreting  the  results  of  these  graphs. 

Let  the  total  run  length  of  N  cylces  be  divided  into  m  replications  of  n 
cycles  per  replication  (or  section).  Let  the  subsection  lengths  within  each  section  be 


{n^,  k  »  1, ...,K}  and  let  be  the  number  of  regenerative  estimates,  AeCj.n^),  formed 

from  n^  cycles,  that  may  be  obtained  from  a  total  of  n  cycles;  then  m^  *  [n/n^] . 

In  the  Retrenched  Graph  (see  e.g.  Figure  3C) ,  the  x-axis  runs  from  0  to  n  cycles. 
Additional  positions  along  the  axis  are  provided  for 

(i)  the  single  *e(mn) ,  the  regenerative  estimate  from  the  whole  sample,  marked  by  * 
and  labelled  REG  on  the  right-hand  side  of  the  graph; 

(ii)  the  single  aA<XAZ(m,n)  estimate,  labelled  ARARE  and  marked  by  o,  with  a 
horizontal  line  through  it  and 

(iii)  the  box  plot  at  position  ARARE  of  the  A.a/Le  estimates. 

Let  OA.e(mm^,njt)  be  the  average  of  the  o'ie(mjt>n^)  's,  i.e.  OAe(mm^,n^)  * 

0/te(j  .m^.n^/m  where  OAe(j  is  the  realization  of  0^(111^,11^)  on  the  jth 

replication.  The  set  of  points  {n^,  ote(mmjt,n^)) ,  k  ■  1,...,K)  are  denoted  by  stars; 
these  are  unbiased  point  estimates  of  (E(A.e(n^)),  k  ■  1,...,K).  The  average  bias  curve, 
i.e.,  the  set  of  points  {(nk>  J*mQ  F^’1),  k  *  1,...,K}  where  d  is  the  degree  of 
fit  and  8^  is  the  average  of  the  m  i.i.d.  estimates  of  the  coefficient  of  n  *  in 
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the  hies  expansion,  is  plotted  as  a  connected  curve.  The  horizontal  line  at  height 
Ohahl( m,n)  represents  the  asymptote  of  the  average  regression  curve.  A  box  plot  of 

either  {aA.e(j,mk,nk) ,  j  -  1 . m)  (as  in  Figure  4G)  or  {^(j.n^),  j  -  l,...,™^} 

(as  in  Figure  4C)  is  displayed  at  x-axis  position  for  each  k.  If  less  than 

nine  points  are  available  for  a  box  plot,  the  individual  values  are  plotted  as  magnitude- 
ordered  circles  rather  than  displaying  a  box  plot  for  so  few  points;  it  is  simple  to 
see  the  spread  and  location  of  the  estimates  this  way. 

A  count  of  the  number  of  hahe( j,n)  estimates  which  are  less  than  Ae(nm)  is 
given:  this  gives  a  rough  sign-test  for  lack  of  bias  in  the  {hxthe(  j  ,n)i ,  with 
Ae(mn)  being  used  as  the  true  value.  The  values  of  aha) lC(m,n),  t>KOKl  (m,n),  Ae(nm) 

A 

and  o(nm),  the  usual  regenerative  standard  deviation  estimate  based  on  all  nm  cycles, 
are  printed  below  the  x-axis.  In  the  BASIC  GRAPH  when  m  ■  1  the  ARARE  standard 
deviation  is  set  equal  to  zero;  the  intent  is  to  eventually  use  here  a  regression  based 
estimate  for  the  standard  deviation  of  Ae(n) : 

With  regard  to  Interpretation,  the  ideal  is  a  Retrenched  Graph  for  which 
(1)  the  average  bias  curve  is  smooth  and,  beyond  [n/2],  close  to  its  asymptote  and 

relatively  flat  ; 

(ii)  the  estimates  he( nm),  ahahe.( m,n)  and  ahz( m,n)  do  not  differ  substantially* 

(iii)  the  distributions  of  {Ae(j,n),  j  ■  l,...,m)  and  {AO»te(j  ,n) ,  j  -  l,...,m) 

are  symmetric  about  Ae(nm)  and  their  respective  means,  and 

(iv)  ihahe.( m,n)  and  o(nm)  are  approximately  equal. 

Figure  4G  is  just  such  a  graph.  When  n  is  small  enough  that  Ae(n)  is  still  biased 
and  hahe.( n)  needs  to  be  used  as  a  less  biased  estimated,  this  bias  shift  will  show 
clearly  in  the  relative  shift  of  the  box  plots  of  the  Ae(j,n)'s  and  the  hahe.(i »n) 's 
as  in  Figure  4C. 

Properties  (1)  and  (ii)  listed  above  and  illustrated  in  Figure  3B  are  the  ideals 
in  a  Basic  Graph.  Because  of  the  high  correlation  between  a*e(2,[n/2])  and  hz( n), 
a  flat  bias  curve  beyond  [n/2]  indicates  that  the  difference  between  E(4£(n))  and 


l°/2])  Is  small.  This  Implies  that  the  bias  In  -te.(n)  is  also  small.  The 
Basic  Graph  is  therefore  very  useful  in  determining  the  section  length,  n,  to  be 
replicated  even  though  the  variability  in  the  estimates  may  still  be  substantial. 

Since  a  Basic  Graph  consists  of  a  single  replication,  only  limited  distributional 
information  can  be  learned  from  it.  However,  the  distributions  of  {^(j.n^), 
j  ■  l,...,m^}  will  give  the  user  a  general  idea  as  to  the  variability  in  the  data. 

Further  interpretation  of  these  graphs  is  given  in  the  protocol  of  Section  5 
and  the  examples  of  Section  6. 

4.3.  Normal  Plots  and  the  Running  Normal  Plot 

To  detect  a  normal  distribution,  a  more  detailed  analysis  of  the  data  than 
is  possible  with  box  plots  is  accomplished  with  normal  plots.  A  number  of  statistical 
tests  accompanying  the  normal  plots  may  be  used  to  formally  test  the  hypothesis  of 
normality.  Normal  plots  [see  e.g.  Wilk  and  Gnanadesikan  (1968)]  are  formed  as  follows. 

Let  X  *  {X^,  X^,  ...  ,  X^}  be  an  i.i.d.  sequence  of  random  variables  with 
continuous  distribution  function  F(x),  -»  <  x  <  “>.  Let  *(x) ,  -*  <  x  <  °°,  be  the 


standard  normal  distribution  function  and  let  and  F-1(*)  be  the  inverse 


distribution  functions  of  $  and  F  respectively.  Let  X^  <  X^  <  <  X 


(n) 


-1, 


be  the  order  statistics  of  X.  An  estimate  of  F  (k/n+1)  is  X^j  and  f°r  any 


.-1. 


constants  a,  b  (b  >  0)  a  plot  of  the  pairs  of  points  {($  (k/n+1),  (X^-a)/b), 


k  ■  1,...  ,  n)  will,  for  large  n,  result  in  a  straight  line  if  X  is  a  set  of 


normally  distributed  random  variables  (i.e.  F  *  4») .  We  standardize  the  plot  by 

2  «2_rn  =.2 


setting  a  -  X  »  ££«i  *yfn  ant*  ^  =^k*l  ^^k~^  /(n-1).  Departure  from  linearity 

in  this  normal  plot  indicates  non-normality. 

Formal  tests  of  normality  that  we  use  include  the  Shapiro-Wilk  statistic 
(1965)  which  was  also  used  by  Fishman  (1977) .  For  ease  of  computation 
we  apply  the  test  to  at  most  50  of  the  points.  The  significance  level  of  the  test 
was  computed  using  the  approximation  in  Shapiro  and  Wilk  [1968].  Additional  test  statis¬ 


tics  are  the  coefficients  of  skewness  and  kurtosis  (Snedecor  and  Cochran  (1967),  p.  86) 
3  ~J  '  '  2'  - -vnJ, 


estimated  by  -  3,  respectively, where  (X^-X^/n. 
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If  the  X^'s  are  normally  distributed,  then  for  large  n  these  moment  statistics  are 
approximately  normally  distributed  with  means  zero  and  standard  deviations  /6/n 
and  v^24/n  respectively.  A  coefficient  of  skewness  or  kurtosis  that  differs  from 
zero  by  more  than  about  two  times  its  standard  deviation  indicates  significant 
departure  from  normality.  Older  tests  of  normality  were  based  on  these  sample 
coefficients;  recently  Interest  has  renewed  is  using  the  pair  of  statistics  as 
a  test  for  normality  (see  Shapiro  and  Wilk  (1965)).  One  can  also  use  Kolmogorov- 
Smirnov  type  statistics  suitably  adapted  for  the  fact  that  the  mean  and  standard 
deviations  are  estimated  from  data  (Lllllefors  (1967)). 

Figures  3E,  3H,  3J,  4F,  and  4H  are  examples  of  a  sequence  of  normal  plots 
of  simulation  generated  estimates;  we  call  this  a  running  normal  plot.  The  goal 
of  these  graphs  is  to  monitor  the  convergence  to  normality  of  the  regenerative 
estimates  Ae(j,n^)  and  to  determine  at  what  point  normality  may  safely  be  assumed 
to  exist.  In  addition  to  the  features  previously  mentioned,  a  line  through  (0,0) 
with  slope  one  is  drawn  for  each  normal  plot  in  the  sequence.  For  any  subsection 
length,  n^,  either  {Ae(j,n^),  j  *  l,...,m^m}  (as  in  Figure  4F  and  labeled  REG)  or 
(aA.e(j,n^),  j  -  l,...,m}  (as  in  Figure  4H  and  labeled  AREG)  may  be  plotted  in 
addition  to  a  plot  of  {AOAe( j,n)  j  *  l,...,m). 

Placing  multiple  normal  plots  on  a  single  running  normal  graph  is  very  useful 
in  identifying  trends  in  convergence.  For  example,  in  Figure  4F  the  plots  of 
{*e(j,500)},  (*e(j , 1000) }  and  {ACW.e( j , 1000) }  when  examined  individually  do  not 
indicate  significant  departures  from  normality.  However,  when  placed  in  sequence 
along  with  the  plot  of  {Ae(j,250)>,  they  are  seen  to  exhibit  a  pattern  of  skewness 
and  nonnormality  similar  to  that  of  {*e(j,250)}  (although  to  a  lesser  extent).  Since 
normal  plots  and  tests  are  difficult  to  Interpret  for  small  data  sets, the  comparison 
between  normal  plots  at  small  and  large  values  of  n  is  helpful. 


I 


5.  A  PROTOCOL  FOR  SEQUENTIAL  APPLICATION  OF  THE  GRAPHICS 
The  following  protocol  outlines  the  sequential  use  of  the  regression  adjusted 
regenerative  estimates  and  graphics  in  analyzing  the  output  of  a  simulation.  The 
objective  is  to  determine  the  total  run  length  needed  to  achieve  a  given  precision  in 
the  simulation  and  to  obtain  a  sectioning  of  that  run  length  into  m  replications  so 
that  the  resulting  4MA.&  estimates  are  approximately  unbiased  and  normally  distributed. 

The  protocol  consists  of  two  stages.  The  first  stage  of  the  sequential  protocol  pro¬ 
vides  a  preliminary  estimate  of  the  section  length  and  the  degree  of  the  fitted  regression 
curve  which  is  needed  to  eliminate  the  concern  over  bias,  i.e.  to  obtain  an  approxi¬ 
mately  bias-free  flOAt  estimate.  It  also  provides  an  initial  estimate  of  the  simula¬ 
tion's  variability.  The  second  stage  of  the  protocol  increases  the  run  length  until 
the  desired  accuracy  is  achieved.  The  degree  and  section  length  are  adjusted  along 
the  way  to  maintain  low  bias  and  Increase  distributional  symmetry.  When  the  final 
accuracy  is  achieved,  further  adjustments  are  made  to  gain  normality  in  the  estimates 
from  each  section. 

We  suggest  that  between  10  and  20  replications  be  available  at  the  end  of  the 
simulation,  since  an  adequate  estimate  of  the  variance  requires  at  least  10  observa¬ 
tions  (replications).  Twenty  observations  should  be  sufficient  for  most  comparison 
purposes  (Mosteller  and  Tukey,  1977,  Ch.  7)  and  keeping  more  than  20  sections  will 
only  slow  convergence  to  normality.  Thus,  for  example,  if  the  plots  showed  that  the 
required  precision  was  obtained  with  40  sections, each  of  length  100  cycles^  and  that 
the  AO/te(j ,100) 's  were  approximately  normal,  we  would  still  recommend  refiguring  the 
simulation  to  20  sections,  each  of  length  200  cycles. 

The  protocol  is  intended  only  as  a  guide  and  a  user  need  not  adhere  to  it 
strictly.  The  whole  point  of  the  graphics  is  to  assist  users  in  exercising  their 
own  judgment.  , 

Throughout  the  protocol  we  use  N  to  represent  the  total  simulation  run 
length  (in  cycles).  The  section  length,  or  number  of  cycles  per  replication,  is 
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denoted  by  n  and  m  Is  the  number  of  sections,  or  replications;  thus  M  ■  [N/n]. 

The  subsection  lengths,  picked  on  the  basis  of  empirical  tests  such  as  those  detailed 
in  Section  7,  are  defined  by  n^  »  [n/k]  for  k  -  1,  2,  3,  4,  7,  10. 

I.  Stage  I-Pilot  run 

1.  Set  Nq  *  initial  run  length. 

2.  Simulate  a  total  of  Nq  cycles. 

a.  Set  m  *  number  of  replications  *  1 

b.  Set  n  *  section  length  *  number  of  cycles /replications  =  Nq 

c.  Set  subsection  lengths  n^  *  [n/k],  k  *  1,  2,  3,  4,  5,  7,  10. 

3.  Form  Basic  Graphs  for  degrees  d  *  1,2,3. 

a.  Pick  smallest  d  so  that  regression  curve  is  smooth  and  fits  data  well.  If  no 
such  d  can  be  determined  increase  Nq  and  go  back  to  1.2. 

b.  Find  n*,  the  smallest  n^  so  that  the  regression  curve  is  nearly  flat  beyond 
n*  and  there  is  little  difference  between  OAZ([ n/n*],n*),  -'UM.e(n)  and 
anedn/n'  ],n')  for  n'  >  n*.  If  no  such  n*  can  be  determined  increase  NQ 
and  go  back  to  1.2. 

c.  Baaed  on  the  spread  of  points  in  the  Basic  Graph  and  the  regenerative 
standard  deviation  estimate,  o(NQ)  project  a  new  estimated  run  length 
N  to  yield  desired  accuracy. 

d.  Set  n  =  n*  and  go  to  Stage  II. 

(Note  that  one  may  conclude  at  this  point  that  the  regenerative  cycles  are  so  long 
that  for  a  given  precision  one  must  accept  a  biased  and/or  non-normal  regenerative 
estimate,  although  it  would  be  better  to  accept  the  less  biased  but  still  non¬ 
normal  regression  adjusted  regenerative  estimate.  Note  too  in  the  case  where  the 
variability  of  ^(Nq)  and  •'UM.e(NQ)  is  high,  the  "power"  of  the  procedures  is 
low  and  the  decisions  made  here  may  have  to  be  modified  in  Stage  II.) 


II.  Stage  II — Retrenchment  and  Fine  Tuning 

1.  Simulate  a  total  of  N  cycles. 

2.  Set  m  ■  [N/n] ,  where  n  was  determined  in  Stage  I. 

3.  Form  Retrenched  Graphs  for  degrees  d  -  1.2,3. 

a.  Pick  smallest  d  giving  good  fit  to  {a/Le(mmk,n^)} . 

b.  If  m  >  20  and  either 

i.  aJ ie(mn)  or  m.n)  show  substantial  bias  as  indicated  by  the 

average  regression  curve  being  far  from  its  asymptote,  or  m.n) 

or  O'te(m.n)  being  significantly  different  from  A&(N),  or 
ii.  Box  plots  of  {/ie(j,n)}  or  {/uM.£(j,n)>,  j  «  l,...,m,  indicate  asymmetry 

then  double  n  and  go  back  to  II. 2. 

c.  Otherwise  if  4/m/i£(m,n)  and  o(N)  are  sufficiently  small  to  yield  the  desired 
accuracy  go  on  to  II. 4;  otherwise  project  a  new  run  length  N  to  yield 
desired  accuracy  and  go  back  to  II.l. 

4.  If  m  >  20  set  m  ■  20;  if  m  <  10  set  m  *  10.  Set  n  *  [N/m] . 

a.  Form  Retrenched  Graph  and  adjust  the  degree  as  in  II. 3. a. 

b.  Form  normal  plots  of  {AiW.e(j  ,n)} ;  {A-etj.n^}  and  {a*e(j  .m^.n^)}  for 
k  -  1,2,4. 

c.  If  normality  is  indicated  then  stop. 

d.  If  normality  is  not  indicated,  then 

1.  if  m*  10  increase  N  if  possible  and  go  to  II.l  ; 
ii.  otherwise  reduce  m,  "keeping  m  >  10,  and  go  to  II. 4. 

Note  that  in  this  Stage  we  are  trying  (1)  to  fix  more  accurately  the 

sample  size  needed  to  attain  the  desired  precision  in  the  final  estimate; 

(ii)  to  section  so  as  to  get  enough  degrees  of  freedom(m-l)  in  the  variance  estimate, 

2 

4 AjCUU.  (m.n),  so  as  not  to  sacrifice  precision  in  this  estimate,  and  (ill)  to  keep  n 
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as  large  as  possible  so  as  to  ensure  that  the  AM£.(J  ,n) ' s  are  approximately 
normally  distributed.  The  decision  may  .  however,  be  that  these  conflicting  objec¬ 
tives  cannot  be  met  and  a  specific  compromise  has  been  adopted .  This  case  will  be 
illustrated  in  the  third  example  in  the  next  section. 

6.  SOME  EXAMPLE  QUEUEING  SYSTEM  SIMULATIONS 
In  this  section  we  give  examples  of  the  use  of  the  protocol  and  the  asso¬ 
ciated  graphics  and  estimates  in  the  analysis  of  the  output  of  some  simple  queueing 
systems.  The  systems  are  two  single  server  queues  with  high  traffic  intensities  and 
a  closed  queueing  network  model  of  a  computer  system.  The  single  server  queues, 
the  M/M/1  and  M/G/l  queues  with  traffic  intensity  0.90  [see  e.g.  Kleinrock,  1975], 
though  simple  in  structure,  exhibit  substantial  bias,  produce  highly  skewed  output 
and  require  very  long  run  lengths  before  stabilizing.  They  therefore  provide 
challenging  tests  of  the  methodology. 

The  closed  queueing  network  on  the  other  hand,  requires  relatively  few 
cycles  before  stabilizing.  This  example  shows  that  the  graphics  can  be  of  use  even 
when  few  cycles  are  needed.  However,  when  the  number  of  cycles  is  so  small  as 
to  preclude  use  of  the  Retrenched  Graph,  the  assumptions  of  low  bias  and  normality 
cannot  be  verified  and  caution  should  be  exercised.  In  fact  the  conclusion  may  be 
the  precautionary  one  that  the  regenerative  estimate  is  very  likely  biased  and  non- 
normal,  the  non-normality  applying  also  to  the  regression  adjusted  regenerative 
estimate. 

a.  The  M/G/l  Queue  (Examples  1  and  2) 

Let  {Wn,  n  >  0}  be  the  waiting  time  of  the  nth  customer  in  a  single 

server  queue  and  let  {A  ,  n  >  1}  and  (S  ,  n  >  0}  be  the  i.i.d.  sequences  of 

n  —  n  — 

interarrival  and  services  times  respectively.  The  waiting  time  sequence  is  defined 

by  Wrt  -  0  and  W  , ,  ■  (W„  +  S  -  A  )+  for  n  >  0  where  x+  »  max(x,0).  Assume 

J  0  n+1  n  n  n+1  — 
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Aq  Is  exponentially  distributed  with  mean  1/X  and  let  SQ  be  hyperexponentially 

distributed;  P(SQ  >  X)  -  l-G(x)  -  p  expC-pjx)  +  (1-p)  exp(-v2x).  Then 

E(Sb)  -  1/p  -  p/Vl  +  (1-p)  A*2*  Variance (Sn)  -  o*  -  2{ph\  +  (l-p)/p2>  '  1/p2 . 

and  the  coefficient  of  variation  of  S  is  C  *  pa  .  If  the  traffic  intensity 

n  s  s 

is  p  ■  X/ji  <  1  then  regenerations  occur  whenever  Wq  •  0  and  the  waiting  time 

2 

has  a  stationary  distribution  with  mean  E(W)  «  p(l  +  C  )/{2p(l-p)}  (see  Kleinrock  (1975)). 

s 

If  p  -  1  then  the  queue  is  known  as  the  M/M/1  queue.  For  the  M/M/1  queue  of 
this  section  (Example  1)  X  *  0.9,  p  ■  1.0,  p  ■  0.9  and  E(W)  «  9.00,  whereas  for 
the  M/G/ 1  queue  (Example  2)  X  *  0.9,  p^  *  0.5,  p2  *  2.0,  p  *  0.33,  p  *  0.9  and 
E(W)  -  13.42. 


b.  Closed  Queueing  Network  (Example  3) . 

Consider  the  queueing  network  pictured  in  Figure  2.  This  network  is  a  model 
of  a  timesharing  computer  system.  There  are  R  service  centers  with  a  fixed  number, 
T,  of  customers  in  the  network.  Service  Center  1  consists  of  T  terminals.  The 
terminal  users  submit  jobs  to  a  computer  system  consisting  of  a  CPU  (Service  Center  2) 
and  a  number  of  peripheral  input-output  devices  (Service  Centers  3,  ...»  R) .  The 
CPU  operates  under  the  processor  sharing  discipline  while  the  peripherial  devices 
are  each  first-come,  first-served  single  server  queues.  Let  p^  be  the  mean  of 
the  exponentially  distributed  service  times  at  service  center  i.  Routing  through 
the  network  is  Markovian  and  the  routing  probabilities  are  given  in  Figure  2 . 

By  assuming  that  ail  service  and  routing  mechanisms  are  mutually  Independent  then 
<£  *  (Q(t),  t  ^  0}  is  a  continuous  time  Markov  chain  where  <£(t)  ■  (Q^(t) , . . .  »Qjj(t)) 
and  Q^(t)  is  the  number  of  customers  in  service  center  i  at  time  t.  The  assump¬ 
tions  on  the  service  distributions  and  disciplines  and  on  the  routing  are  such  that 
the  equilibrium  distribution  of  the  network  exists  and  has  a  product  form  (Baskett, 
et  al.  (1975)).  Define  a  response  time  to  be  the  time  from  when  a  customer  leaves 
service  center  1  until  that  customer  next  returns  to  service  center  1  and  let 

n 

be  the  just  completed  response  time  of  the  nth  customer  to  arrive  at  service  center  1. 


i 

1 

,i 

1 

J 

i 


4 
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Then  W  -  (W^,  n  _>  0}  is  regenerative  with  regenerations  occurring  at  points  n 
such  that  the  nth  customer  arrives  at  service  center  1  leaving  centers  2,  ...  ,  R, 
empty.  Again  we  shall  be  interested  in  the  expected  stationary  response  time 
E(W)  which  is  known  to  be  finite  (Lavenberg  and  Sauer  (1977)).  Let  be 

the  steady  state  utilization  of  service  center  i.  The  particular  parameters  chosen 
for  this  model  are  listed  in  Figure  2  and  yield  ?2  m  0.894,  *  0.268  for  i  >.  3 

and  E(W)  -  8.65. 

The  application  of  the  methodology  to  these  three  examples  is  given  in  the 
Figures;  each  Figure  has  a  Caption  and  an  Interpretation  which,  when  read  sequen¬ 
tially  describe  the  evolution  of  the  simulation. 


Closed  Queueing  Network 


1 


T  -  10 

y'1  -  10.0 

p  -  0.850 
P2  -  0.895 


System  Parameters 


R  = 

6 

E(W)  *  8.650 

-1 

y2  = 

0.250 

u"1  -  0.300  for  i  >  3 

Pi  " 

0.250 

for 

i  >  3 

P^  * 

0.268 

for 

i  >  3 

FIGURE  2.  Closed  queueing  network  which  models  a  timesharing  computer 
system. 
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BASIC  GRAPH 
M/M/1  p=0 . 9 

1  REPLICATIONS  OF  500  CYCLES 


ARARE  ESTIMATE  :  10.28 

ARARE  STD  DEV  :  .00 

REG  ESTIMATE  :  10.13 

REG  STD  DEV  :  2.08 

DEGREE  OF  FIT  :  3 


FIGURE  3 A.  Example  1.  Stage  I,  steps  1  and  2.  No  set  at  500  cycles.  A 
decision  has  been  made  on  the  basis  of  similar  graphs  that 
d  »  3  is  probably  needed. 

Interpretation.  Since  the  precision  (2o/u)  is  about  40Z,  and  the 
regression  curve  is  still  coming  up  to  the  asymptote.  No  •  500 
may  be  too  small.  Therefore  return  to  Stage  1.1  and  set 
Nq  -  1000. 
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BASIC  GRAPH 
M/M/1  p=0 . 9 

1  REPLICATIONS  OF  1000  CYCLES 


ARARE  ESTIMATE  :  10.74 

ARARE  STD  DEV  :  .00 

REG  ESTIMATE  :  10.71 

REG  STD  DEV  :  2.03 

DEGREE  OF  FIT  :  3 


FIGURE  3B.  Example  1.  Stage  I,  steps  1  and  2  with  Nq  Increased  to  1000. 

Interpretation .  Since  the  regression  curve  beyond  500 
and  *e(1000)  and  1000)  are  approximately  equal,  simulation 

out  to  n  ■  1000  is  likely  to  be  adequate.  Since  it  is  estimated 
that  about  16  times  as  many  cycles  are  needed  to  cut  the  standard 
deviation  down  by  a  quarter,  an  exploratory  excursion  to 
Stage  II  with  N  ■  4000  is  taken. 
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RETRENCHED  GRAPH 
M/M/1  p= 0.9 

8  REPLICATIONS  OF  500  CYCLES 
BOXPLOTS  OF  AREG  ESTIMATES 


FIGURE  3C.  Example  1.  Stage  II,  step  3a.  N  »  4000  and  d  cut  back  to  2. 

Also  n  is  cut  back  to  500  to  determine  whether, with  this 
greater  precision,  a  firm  decision  can  be  made  on  section  length  n. 
Note  that  the  box  plots  are  of  the  aAe(j,nk)  which  show  less  skew¬ 
ness  then  the  AC(j,nk)  in  Figure  3b. 

Interpretation.  The  4/Ut*e(8,500)  of  0.78  is  lower  than  the 
regenerative  s.d.  estimate  of  0.91.  Clearly  about  two  to  three 
times  as  many  cycles  are  needed.  There  is  also  still  bias,  as 
illustrated  by  the  upward  displacement  of  the  data  points  at  ARARE 
compared  to  the  data  points  at  n  ■  500.  This  shift  is  not 
statistically  significant  but  is  reinforced  by  the  fact  that  the 
regression  curve  is  increasing.  The  degree  d  ■  2  is  now  seen 
to  give  a  good  fit.  For  illustrative  purposes  the  next  step  is 
to  repeat  this  Retrenched  Graph  with  twice  the  number  of  cycles. 
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RETRENCHED  GRAPH 
M/M/1  p=Q . 9 

16  REPLICATIONS  OF  500  CYCLES 


ARARE  ESTIMATE  :  9.76 

ARARE  STD  DEV  :  .54 

REG  ESTIMATE  :  9.72 

REG  STD  DEV  :  .65 

DEGREE  OF  FIT  :  2 


FIGURE  3D.  Example  1.  Stage  II,  step  3  repeated  with  N  -  8000  for  illustration. 

Interpretation.  The  presence  of  slight  bias  in  A£(j,500)  is 
confirmed.  It  is  also  clear  from  the  standard  deviation  estimates  that 
more  cycles  are  needed  to  reduce  the  standard  deviation  to  less  than 
0.5.  At  this  point  a  running  normal  plot  is  helpful  and  is  given  in 
the  next  figure. 


M/M/1  p=0 . 9  NORMAL  PLOTS 
16  REPLICATIONS  OF  500  CYCLES 


REGC  125] 
SHAPIRO-WILK  . 866(  .000) 
SKEW  1 . 143 (  .306) 
KURTOSIS  . 693(  .612) 


REGC  250]  REGC  500] 

. 955(  .234)  . 962 (  .673) 

.362 (  .433)  - . 216  (  .612) 

-.751 (  .866)  -.708(1.225) 


RARE 

. 974(  .873)  (SIG  LEV) 
-  . 075 (  .612)  (SD  SKEW) 
-.634(1.225)  (SD  KURT) 


FIGURE  3E.  Example  1.  Supplement  ot  Figure  3D.  N  ■  8000  cycles. 

Interpretation.  Note  that  while  the  ^tE(jt125),s  are  clearly  non¬ 
normal,  there  Is  no  Indication,  graphical  or  from  the  test  statistics, 
of  departure  from  normality  either  In  -*te(j,500)  or  fiaAt( j,500).  The 
number  of  cycles  will  now  be  doubled  to  reach  the  desired  precision; 
but  first  the  need  to  go  to  n  ■  1000  for  reasons  of  bias  Is  explored. 
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RETRENCHED  GRAPH 
M/M/1  p= 0.9 

8  REPLICATIONS  OF  1000  CYCLES 


ARARE  ESTIMATE  :  9.70 

ARARE  STD  DEV  :  .63 

REG  ESTIMATE  :  9.72 

REG  STD  DEV  :  .65 

DEGREE  OF  FIT  :  2 


FIGURE  3F.  Example  1.  N  -  8000  cycles.  Illustrative  retrenched  graph  giving 
confirmation  of  the  basic  graph  in  Figure  3b. 

Interpretation.  The  bias,  if  any,  at  n  ■  1000  is  small,  and 
would  be  taken  care  of  by  using  the  ,1000) '  s  instead  of  the 

>te(j ,1000) 's.  Notice  the  evolution  of  the  extreme,  low  outlier. 
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RETRENCHED  GRAPH 
M/M/1  p=0 . 9 

32  REPLICATIONS  OF  500  CYCLES 
BOXPLOTS  OF  AREG  ESTIMATES 


ARARE  ESTIMATE  :  9.11 

ARARE  STD  DEV  :  .35 

REG  ESTIMATE  :  9.06 

REG  STD  DEV  :  .42 

DEGREE  OF  FIT  :  2 


FIGURE  3C.  Example  1.  Stage  II,  step  3  with  N  ■  1600  cycles. 

Interpretation.  The  two  estimated  standard  deviations  are  close 
and  indicate  that  at  least  10%  precision  has  been  attained.  We 
are  at  Stage  II,  Step  3C  here  and  ready  to  go  to  Stage  II,  Step  4. 
Since  m  >  20  and  there  is  some  bias,  n  should  at  least  be  doubled. 
First  we  look  at  the  normality  in  the  next  figure. 
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M/M/l  p= 0.9 
32  REPLICATIONS 


NORMAL  PLOTS 
OF  500  CYCLES 


o 


REGC 

1253 

REGC 

250] 

REGC 

500] 

RARE 

SHAPIRO-WILK 

.866( 

.000) 

.  942  ( 

.029) 

.961  ( 

.339) 

.  965  ( 

.428) 

(SIG  LEV) 

SKEW 

1 . 124 ( 

.217) 

.592( 

.306) 

.  265  ( 

.433) 

.  370  ( 

.433) 

(SD 

SKEW) 

KURTOSIS 

1 . 057 ( 

.433) 

'T 

O 

■'T 

i 

.612) 

- . 829 ( 

.866) 

- . 570 (  .866) 

(SD 

KURT) 

FIGURE  3H.  Example  1.  N  *  16000  cycles. 

Interpretation.  Another  look  at  normality  with  twice  the  amount  of  data 
as  was  available  in  Figure  3D.  Again  there  is  no  formal  indication  of  depar¬ 
ture  from  normality  in  the  AOAe(j ,500) 's,  but  it  is  to  be  recalled  that  with 
only  32  data  points,  tests  for  normality  have  relatively  low  power.  A  visual 
comparison  of  the  AOA.e(j ,500) ’s  to  the  Ae(j ,125) ’s  and  Ae(j ,250) 's  indicates,  k 
however,  the  possibility  of  nonnormality.  Since  m  *  32  >  20,  n  is  now  doubled. 
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RETRENCHED  GRAPH 
M/M/1  p= 0.9 

16  REPLICATIONS  OF  1000  CYCLES 


100  332  1000  ARARE  REG 


ARARE  ESTIMATE  : 

9.06 

ARARE  STD  DEV  : 

.37 

REG  ESTIMATE  : 

9.06 

REG  STD  DEV  : 

.42 

DEGREE  OF  FIT  : 

2 

FIGURE  31.  Example  1.  N  *  16,000  cycles.  Final  steps  of  Stage  11  of  the 

protocol.  Note  the  similarity  between  this  Figure  and  Figure  3F 
illustrating  the  stability  reached  with  this  sample  size. 

Interpretation.  The  bodies  of  the  box  plots  of  the  Aa/te(j  ,1000)  's 
and  the  AC(J ,1000) *s  (for  J  ■  1,...,16)  indicate  that  there  is  very 
little  to  choose  between  them;  their  spreads  are  similar.  Also  there 
are  no  extreme  values  or  outliers.  The  estimated  standard  deviation 
of  0.37,  MOAe(16,000)  indices  that  better  than  10Z  precision  has 
been  achieved.  A  retrenched  graph  with  degree  d  -  3  is  indis¬ 
tinguishable  from  this  graph. 
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M/M/1  p=0 . 9  NORMAL  PLOTS 
16  REPLICATIONS  OF  1000  CYCLES 


REGC 

250] 

REGC 

5003 

REGC  10003 

RARE 

SHAPIRO-WILK 

.941  ( 

.026) 

.961  ( 

.339) 

. 992(  .999) 

. 988 (  .993) 

(SIG  LEV 

SKEW 

.592( 

.306) 

.  265  ( 

.433) 

- . 066 (  .612) 

- . 007 (  .612) 

(SD  SKEW 

KURTOSIS 

i 

o 

.612) 

- . 829 ( 

.866) 

-.383(1.225) 

-.319(1.225) 

(SD  KURT 

FIGURE  3J.  Example  1.  N  *  16,000  cycles.  Final  look  at  normality  Stage  II,  Step  4C. 

Interpretation .  For  illustration  the  *e(j .n^) 's  are  shown,  rather 
than  the  OA.fe(j,njc).  There  is  no  indication  of  departure  from  normality 
in  the  *O*e(j,1000)  's  or  the  A.e(j  ,1000) 's. 


BASIC  GRAPH 
M/G/l  p=0 . 9 

1  REPLICATIONS  OF  500  CYCLES 


50  250  500  ARARE  REG 


ARARE  ESTIMATE  :  11.06 

ARARE  STD  DEV  :  .00 

REG  ESTIMATE  :  10.37 

REG  STD  DEV  :  2.76 

DEGREE  OF  FIT  :  3 


FIGURE  4A.  Example  2.  Stage  1,  steps  1  and  2.  Nq  set  at  500  cycles. 

A  decision  has  been  made  on  the  basis  of  similar  graphs  that 
d  -  3  is  probably  needed. 

Interpretation.  Because  of  the  skewed  service  time,  there 
is  more  variability  in  this  queue  than  in  the  M/M/1  queue.  The 
regression  curve  is  definitely  still  rising  so  that  n  ■  500 
may  be  too  small.  Also  the  precision  is  about  50%,  so  that  about 
25  times  more  cycles  are  needed.  An  exploratory  Retrenched  Graph 
is  needed. 
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RETRENCHED  GRAPH 
M/G/l  p=0 . 9 

4  REPLICATIONS  OF  500  CYCLES 
BOXPLOTS  OF  REG  ESTIMATES 


50  250  500  ARARE  REG 


ARARE  ESTIMATE  :  11.72 

ARARE  STD  DEV  :  .74 

REG  ESTIMATE  :  10.89 

REG  STD  DEV  :  .94 

DEGREE  OF  FIT  :  3 


Figure  4B.  Example  2.  Stage  II,  step  3a.  NQ  -  2000  and  n  -  500. 

Interpretation.  There  Is  clearly  still  bias  in  the 
*e(j ,500) *s  and  a  need  for  more  precision.  Four  times  as 
many  cycles  are  possibly  needed;  this  is  the  "otherwise"  of 
Stage  II,  Step  3c. 


RETRENCHED  GRAPH 
M/G/l  p=0 . 9 

LB  REPLICATIONS  OF  500  CYCLES 
_  BOXPLOTS  OF  REG  ESTIMATES 


•  m  • 


ARARE  ESTIMATE  :  12.77 

ARARE  STD  DEV  :  1.01 

REG  ESTIMATE  :  12.70 

REG  STD  DEV  :  1.04 

DEGREE  OF  FIT  :  3 


Figure  4C.  Example  2.  Stage  II,  Steps  1  through  3  with  NQ  -  8000  and  n  -  500 

Interpretation.  Note  that  the  estimated  standard  deviations 
have  Increased  from  their  values  in  Figure  4B.  This  is  because 


RETRENCHED  GRAPH 
M/G/l  p=0.9 

32  REPLICATIONS  OF  500  CYCLES 
BOXPLOTS  OF  AREG  ESTIMATES 


SO  250  500  ARARE  REG 


ARARE  ESTIMATE  :  13.62 

ARARE  STD  DEV  :  .94 

REG  ESTIMATE  :  13.78 

REG  STD  DEV  :  .87 

DEGREE  OF  FIT  :  3 


FIGURE  4D.  Example  2,  Stage  II,  Steps  1  through  3  with  N.  -16000  and 
n  -  500.  0 

Interpretation.  The  point  estimates  have  risen  in  value 
and  the  estimated  standard  deviations  have  dropped,  but  not 
quite  enough  to  give  10%  precision.  Note  that  the  true 
E(W)  -  13.42.  There  is  clearly  still  some  bias  in  the 

>500) ' s ;  of  the  ,500) 's,  nineteen  are  less  than  or 

equal  to  4e( 16000) . 
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RETRENCHED  GRAPH 
M/G/l  p= 0.9 

20  REPLICATIONS  OF  1000  CYCLES 
BOXPLOTS  OF  AREG  ESTIMATES 


FIGURE  4E.  Example  2.  Stage  II,  steps  1  through  3  with  Nq  increased  to 
20,000,  the  maximum  number  which  computing  constraints  allowed. 
The  degree  of  fit  is  d  ■  2,  sufficient  for  the  purpose,  as 
compared  to  d  ■  3  in  Figure  4D.  The  section  length  n  is 
increased  to  n  ■  1000. 

Interpretation.  There  is  still  indication  of  bias  at 
n  ■  1000,  but  more  particularly  of  outliers  in  the  Kt  (j ,1000) 's 
and  /WA.e(j  ,1000)  '  s  which  need  to  be  looked  at  in  a  running 
normal  plot.  Note  that  it  is  not  quite  possible  to  attain  10% 
accuracy. 
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M/G/l  p=0 . 9  NORMAL  PLOTS 
20  REPLICATIONS  OF  1000  CYCLES 


o 

© 


REGt 

2503 

REGC 

500  3 

REGC  10003 

RARE 

SHAPIRO-WILK 

.890  ( 

.000) 

.  9S8( 

.409) 

.  959  ( 

.516) 

. 943(  .265) 

(SIG  LE\ 

SKEW 

1 . 142 ( 

.274) 

.  487  ( 

.387) 

.606( 

.548) 
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(SD  SKE. 

KURTOSIS 

1 . 154  ( 
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- . 097 ( 

.775) 

-.171(1.095) 

. 056(1. 095) 

(SD  KlilT 

FIGURE  4F.  Example  2.  Running  normal  plot  accompanying  Figure  4G. 

Interpretation.  The  definite  non-normality  in  the 
,250) ' s  is  no  longer  apparent  in  the  values  of  the  test 
statistics  at  n  ■  1000.  However  outliers  are  clearly  apparent 
in  the  *e(J,1000) 's  and  >UW.e(j ,1000) 's,  suggesting,  since 
m  *  20,  that  n  be  doubled.  (Stage  II,  step  4d.) 


RETRENCHED  GRAPH 
M/G/l  p= 0.9 

10  REPLICATIONS  OF  2000  CYCLES 
BOXPLOTS  OF  AREG  ESTIMATES 


ARARE  ESTIMATE  :  13.54 

ARARE  STD  DEV  :  .83 

REG  ESTIMATE  :  13.44 

REG  STD  DEV  :  .75 

DEGREE  OF  FIT  :  2 


FIGURE  4G.  Example  2.  Stage  II,  step  4  with  n  increased  to  2000, 
which  keeps  m  2  10,  as  required.  To  achieve  the  desired 
precision  4000  more  cycles  would  have  been  desirable. 

Interpretation.  The  bias  has  been  reduced  but  more 
importantly  the  heavy  tails  in  the  estimates  have  disappeared 
(as  compared  to  Figure  4E)  at  all  n^'s.  To  confirm  these 
results  a  running  normal  plot  is  created. 


37 


P3 


M/G/l  p=0 . 9  NORMAL  PLOTS 
10  REPLICATIONS  OF  2000  CYCLES 


SHAPIRO-WILK 

SKEW 

KURTOSIS 


AREGC  5001  AREGC 10003  AREGC 2000]  RARE 

. 984 (  .982)  . 971 <  .890)  .9 77 (  .944)  .979(  .952) 
.319(  .775)  . 143(  .775)  . 1 18 (  .775)  .032(  .775) 

- . 671 ( 1 . 549)  - 1 . 090 ( 1 . 549) - 1 . 036 ( 1 . 549) -1.048(1. 549) 


(SIG  LEV) 
(SD  SKEW? 
(SD  KURT) 


FIGURE  4H.  Example  2.  Final  check  for  normality  with  n  ■  2000  cycles. 

Note  that  the  plots  at  500  and  1000  are  of  the  <ViaAZ(j  .n^)  's. 

Interpretation.  The  test  statistics  and  the  normal  plots 
confirm  the  apparent  normality  of  the  estimates  at  n  «  2000, 
although  the  sample  size  m  -  10  is  very  small.  The  estimate 
(VuVLZ(10,  2000)  would  be  preferred  to  O*e(10,2000)  because 
of  bias  considerations. 
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BASIC  GRAPH 
QUEUEING  NETWORK  1 
1  REPLICATIONS  OF  50  CYCLES 


ARARE  ESTIMATE  :  8.43 

ARARE  STD  DEV  :  .00 

REG  ESTIMATE  :  8.49 

REG  STD  DEV  :  .23 

DEGREE  OF  FIT  :  2 


FIGURE  5A.  Example  3.  Stage  I,  steps  1,  2  and  3.  NQ  set  at  50  cycles: 

Degree  of  fit  is  d  ■  2. 

Interpretation.  This  initial  basic  graph  shows  that  the  bias  is 
probably  very  low  at  n  ■  50  and  that  the  accuracy  is  already  less  than 
6%.  Clearly  the  A.e(j,5)'s  are  non-normal  but  there  is  insufficient 
data  to  judge  the  distributions  of  the  other  estimates.  Since  cycles 
are  relatively  expensive,  in  practice  a  decision  to  terminate  the 
situation  at  this  point  would  be  made.  Thus  the  assumptions  of  low 
bias  and  normality  are  not  confirmed.  For  purposes  of  exposition, 
however,  the  50  cycles  are  replicated  25  times  for  the  next  figure. 


RETRENCHED  GRAPH 


QUEUEING  NETWORK  1 
-.25. REPLICATIONS. OF.  50  CYCLES 

BOXPLOTS  OF  REG'  ESTIMATES 


0  I _ LJ I — I - L- 

5  16 


1 


_J _ I - 1 

50  ARARE  REG 


ARARE  ESTIMATE  : 

8.64 

ARARE  STD  DEV  : 

.08 

REG  ESTIMATE 

8.68 

REG  STD  DEV  : 

.06 

DEGREE  OF  FIT  : 

2 

FIGURE  5B.  Example  3.  Stage  II.  N  •  1250. 

Interpretation .  This  retrenched  graph  shows  clearly  that  bias 
is  nonexistent  at  n  -  50.  Note  that  16  of  the  j,50)'s  are 

less  than  4£(1250).  Note  too  the  extreme  negative  skewness  of  the 
AC(j,50)'s;  this  seems  to  have  disappeared  by  n  »  50  both  for  the 
*e(j,50)  's  and  ,50)  1  s.  The  very  flat  bias  curve  and  the 

slight  asymmetry  in  the  distribution  of  the  ,50)  's  suggest 

that  the  degree  of  fit  in  the  regression  could  be  cut  back  to  one. 
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7.  EMPIRICAL  RESULTS 


In  this  section  we  empirically  compare  the  relative  properties  of  the  estimate 
*e(mn)  and  <Mux/te(m,n)  for  the  single  server  queues.  These  studies  give  the  results 
alluded  to  in  Section  3;  namely  that  the  0/t0/te(m,n)  estimate  has  low  bias,  and  variance 
and  mean  squared  error  (mse)  comparable  to  those  of  the  regenerative  estimate. 
Furthermore  confidence  intervals  formed  using  a'ia/ie(m,n)  and  te(m,n)  are  valid, 
i.e.  the  confidence  Intervals  contained  the  true  value  with  the  correct,  prespecified 
probability.  Thus  even  for  large  samples,  a  user  looses  nothing  by  switching  from 
the  regenerative  to  the  CUiaAZ  estimate.  The  substantial  benefit  of  doing  so  is  that 
the  user  may  easily  judge  what  constitutes  a  "large  sample"  and  at  what  point  the 
asymptotic  theory  is  valid.  In  addition  <VuVLe.(m,n)  and  m,n)  were  less 

correlated  with  one  another  than  were  the  regenerative  point  and  standard  deviation 
estimates,  4e(mn)  and  o(nm). 


The  above  mentioned  experiments,  the  results  of  which  are  compiled  in 
Tables  1,  2,  and  3,  were  conducted  as  follows.  The  systems  tested,  as  in  Examples  1 
and  2,  were  highly  congested  M/M/1  and  M/G/l  queues  (p  »  0.9).  For  each  system 
and  several  combinations  of  m  sections,  n  cycles  per  section  and  subsection  lengths 
{n^},  R  independent  simulations  of  m  x  n  cycles  were  performed.  These  queues  were 
simulated  using  the  FORTRAN  language  and  the  random  number  generator  described  in 
Learmonth  and  Lewis  (1973).  Different  seeds  for  the  generator  were  used  for  differ¬ 
ent  systems  and  different  values  of  m  x  n.  On  each  replication  the  following 


quantities  were  computed:  *e(nm),  <&*WAe(m,n),  a/tefmm^.n^) ,  a(mn),  4>uW.e(m,n)  and 
Aa%C(mn^,n^) .  The  entries  in  the  tables  are  best  explained  by  an  example. 

Consider  the  first  row  of  Table  2.  Let  4e^(4000)  and  5^(4000),  for 
1  -  1, . . . ,R,  be  respectively  the  realizations  of  the  regenerative  point  and  standard 
deviation  estimates  based  on  the  4000  cycles  of  the  ith  replication  (m  ■  8,  n*  500). 
The  column  labeled  W  contains  the  sample  average,  say  W^,  of  (*e^(4000), 
i  ■  1,...,R);  ■  ^i«l  **  (^000)/R  »  8.76.  This  is  an  unbiased  estimate  of 
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E(*C(4000)) »  whose  value  Is  unknown.  The  sample  standard  deviation  of  Wfc/1  is 
"  ^i-l{>lei(4000)  "  V2/*®-1)]172  ■  °*09  and  18  listed  in  the  colum 
labeled-  s(W).  'Thus  -a 90%- confidence*  interval  foe-  •  i(*a(4000)->  is  -8.76'+  1.645*0.09- 
»  8.76  +  0.15.  Since  for  this  example  E(W)  ■  9.00  does  not  fall  within  this  con¬ 
fidence  Interval  we  conclude  that  4fc(4000)  is  significantly  biased.  An  estimate 
of  the  mean  squared  error  of  *.£(4000),  MSE^  *  (Ae^(4000)  -  E(W))2/R  *  0.87, 

is  listed  in  the  MSE  column  of  Table  2.  Similarly  estimates  of  the  correlation 
and  Spearman,  or  rank,  correlation  coefficients  (Snedecor  and  Cochran  (1967)) 

A 

between  -42(4000)  and  a (4000),  0.77  and  0.67,  respectively,  are  listed. 


Using  ^(4000)  and  o^(4000)  alleged  90%  confidence  intervals, 

A 

4^(4000)  +  1.645o^(4000) ,  for  E(W)  were  also  formed.  The  fraction  of  these 
confidence  intervals  that  actually  contained  E(W)  is  reported  in  the  column  labeled 
"90%  coverage."  This  fraction,  called  "(90%)  coverage,"  should  be  close  to  0.90 
if  in  fact  valid  confidence  intervals  are  being  formed.  Any  coverage  less  than 
or  equal  to  0.85  in  Table  2  (with  R  »  100)  and  any  coverage  less  than  or  equal 
to  0.865  in  Tables  1  and  3  (with  R  ■  200)  are  significantly  less  than  0.90 
at  the  0.10  level. 


Table  1  focuses  on  the  effects  of  changing  the  set  of  subsection  lengths 
(n^)  on  the  (UUUiZ  estimate.  The  table  reports  results  of  simulations  of  2000 
cycles  of  the  M/M/1  queue  with  p  -  0.9.  The  graphs  of  Section  6  show  this  to  be 
a  relatively  short  run  length.  The  estimates  in  this  table  used  a  polynomial 
of  degree  3  in  the  regression;  degrees  1  and  2  are  omitted  since  they  showed  the 
same  trends  but  with  higher  bias  and  lower  coverage.  Notice  that  the  (VUViZ 
estimates  for  which  the  smallest  value  of  n^  is  n/10  (subsection  sets 


S^,  S^,  S i  and  Sg  in  the  table)  show  less  bias  and  higher  coverage  than  those 
with  values  of  n^  <  n/10  (sets  S^,  and  Sg) .  The  values  of  n^  in 

subsection  sets  Sj  and  S?  are  n/10,  [n/7],  n/5,  n/4,  [n/3],  n/2,  and  n, 
which  is  the  recommended  set  of  points.  Due  to  the  high  correlation  between  the 
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estimates  for  different  values  of  n^,  there  is  little,  if  any,  advantage  in  using 
.  jaqre  values  of  n^  than  are  in  the  recommended  sets  (compare  sets  versus 
and  versus  Sg) .  The  <VUVit  estimate  is  relatively  insensitive  to  the  subsection  '  *  * ' 
lengths  as  long  as  they  are  not  too  small.  However,  in  some  cases  adding  points  that 
are  not  divisors  of  .n  will  increase  the  aJUXAt'a  variance.  The  set  of  points  recom¬ 
mended  above  is  relatively  small  and  has  consistently  yielded  satisfactory  re8ults. 

Tables  2  and  3  report  the  results  of  regenerative  simulations  with  more 
cycles  than  the  simulation  of  Table  1.  The  main  trends  to  notice  in  these  tables 
are  that: 


(i)  the  bias,  standard  deviation,  coverage  and  mse  of  0ta't£.(m,n)  is  comparable 
to  that  of  fil( mn)  ; 

(ii)  a’Ufte(m,n)  shows— lover  bias,  lower  mse  and  higher  (truer)  coverage  than  OAe(m,n), 
(ill)  for  fixed  values  of  m  x  n  it  is  preferable  to  have  small  values  of  m, 
the  number  of  sections,  in  terms  of  bias,  mse,  coverage  and  correlation 
between  point  and  standard  deviation  estimates,  and 

A 

(iv)  the  correlation  between  *e(mn)  and  o(mn)  is  higher  than  that  between 
<VUVilL{ m,n)  and  4A0Ae.(m,n) . 

This  last  property  warrants  further  discussion.  As  gets  large  the  o/te(j ,m^,n^> *s 

converge  to  normal  random  variables  and  since  sufie.(j,n)  is  a  linear  combination  of 
{^^(j.m^.n^) ,  k  -  1,...,K),  it  too  converges  to  normality.  Thus  aAOAe(m,n) 
and  4AAte( m,n)  are  asymptotically  Independent.  The  nonzero  correlations  in 
Tables  2  and  3  indicate  that  for  these  highly  congested  queues  the  convergence 
to  normality  is  very  slow  and  that  the  independence  of  mean  and  variance  is  quite 
sensitive  to  the  normality  assumption.  Of  course  these  correlations  are  seen  in  ■ 

a  J 

Tables  2  and  3  to  be  less  than  the  correlation  between  Ae(mn)  and  o(mn).  Also 
non-parametric  methods  can  be  used  with  the  set  of  m  40Ae(j,n)'s. 
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The  high  correlation  between  Ae(N)  and  a(N),  where  N  *  mn,  is  of  particular. 


concern  since  sequential  rules  for  determining  simulation  run  lengths  typically  rely 

on  the  relative  width  of  a  confidence  interval,  a  multiple  of  o(N)/Ae(N),  as  a 

. .  .  .  •  . .  . . 

stopping  criterion  [Lavenberg  and  Sauer  (1977)].  If  Ae(N)  is  smaller  than 
usual,  then  because  of  the  high  correlation  o(N)  will  also  be  smaller  than  usual. 
This  leads  to  a"low  miss" of  the  true  answer,  i.e.  an  unusually  small  confidence 
Interval  centered  about  an  unusually  small  point  estimate.  Since  sequential  rules 
are  specifically  designed  to  stop  on  small  confidence  intervals,  the  rule  inevitably 
stops  a  large  proportion  of  low  misses.  This  will  not  be  a  problem  with  the  graphical 
methodology  given  here,  since  in  addition  to  the  lower  correlation  between  <MuM.e(m,n) 
and  44aA&(m,n) ,  the  distributional  information  displayed  will  guard  against  the 
skewed  data  which  causes  the  low  miss. 


8.  SUMMARY  AND  FURTHER  WORK 

The  statistical  and  graphical  methodology  given  in  this  paper  has  been  shown 
to  be  very  effective  in  verifying  the  bias  and  normality  properties  of  regenerative 
estimates  and  in  sectioning  up  a  regenerative  simulation  so  that  more  reliable  standard  i 
deviation  estimates  and  confidence  interval  estimates  can  be  obtained  than  with  the  ] 

usual  regenerative  methodology. 

Although  the  design  of  the  regression  for  the  regression-adjusted  regenerative  ] 

j 

estimates  has  been  shown  to  be  robust  in  the  cases  considered,  analytic  results  to  j 

confirm  this  in  a  broader  class  of  simulation  situations  will  be  pursued.  The  j 

regression-based  variance  estimate  for  the  regenerative  estimate  also  needs  to  be  j 

\ 

explored  further.  This  will  be  useful  when  applications  of  the  methodology  to  other  j 

statistical  and  simulation  output  situations  is  considered.  J 

! 

Another  point  to  be  explored  is  the  use  of  robust  regression  techniques, 
as  at  311.  This  could  be  fruitful  because  the  AeCj.n^'s  are  highly  non-normal,  j 

as  can  be  seen  in  Figure  AC  at  n^  •  50.  . 


Finally  we  note  that  the  regression-adjusted  estimation  procedure  and  its 
associated  graphics  can  be  applied  in  many  other  areas  of  statistics »  and  in 
particular  in  the  statistical  analysis  of  simulation  outputs.  Thus  if  one  has  an 

A  A 

estimator  0(n)  ■  0(X. )  of  a  parameter  0  from  an  i.i.d.  sample  X.,,...,X  , 

in  in 

similar  estimates  0(j(n^)»  j  -  1 . m^,  can  be  formed  from  m^  subsections  of 

the  sample  for  k  -  1,...,K  to  give  a  regression-adjusted  estimator,  with  graphics, 
just  as  it  was  done  for  the  regenerative  estimator  A£(n).  This  application, 
and  multivariate  extensions,  will  be  explored  elsewhere. 
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Table  1 

Effects  of  Subsection. lengths  on  A tLtvut  Estimates  of.  E (U)» «  .9 . 00. 
in  the  M/M/1  Queue  with  p  ■  0.9.  Here  R  ■  200. 


Estimate 

Set  of 
Subsection 
Lengths  for 

A naAe. 

Ae<2000) 

<t*ie  (io,  200) 

aie(20,i00) 

aiaie( 10,200) 

sl 

aidie(l0,200) 

S2 

aiaie(10,200) 

S3 

a/ta/te( 10,200) 

S4 

aKne(20,i00) 

S5 

aiaie(20,100) 

S6 

aiaieC 20,100) 

S7 

(futle(20,100) 

S8 

W 


8.725 

8.086 

7.513 

8.466 

8.459 

8.753 

8.673 

8.172 

8.129 

8.466 

8.459 


Degree  of  polynomial  fitted  in  regression  was  3. 

Sj  ■  {i^sn^  "  5k,  k  -  1, . . .  ,40} 

52  "  {nk:nk  *  5,10»  20’  28‘  40>  50»  66»  100»  200} 

53  -  -  20,28,  40,  50,  66,100,200} 

■  {n^n^  *  20  +  (k-l)5,  k  -  1,...,  37} 

55  -  (xyi^  -  4  +  (k-l)3,  k  -  1 . 33} 

56  "  {nk:nk  "  4*  7’  10»  14,  20,  25,  33’  50»  100} 

57  “  {V°k  "  10,14,  20,  25,  33,  50,100} 

So  -  {n^n,.  -  10  +  (k-l)3,  k  -  1,...,  31} 


90% 

Coveraee 

MSE 

.810 

1.552 

.715 

1.847 

.460 

2.948 

.820 

1.707 

.820 

1.424 

.865 

1.661 

.850 

1.659 

.725 

1.864 

.730 

1.691 

.790 

1.485 

.780 

1.592 

Table  2 


Simulation  Results  for  Estimates  of  E(W)  -  9.00 

.  .  ^  .  «  .  •  «  * 

In  M/M/1  Queue  with  p  •  0.9.  Here  R  -  100. 


Estimate 

-  r 

W 

■ 

90% 

Coverage 

MSE 

Correlation 

(W.o(W)) 

Spearman 

Correlation 

(W.o(W)) 

*e(4000) 

8.76 

.09 

.84 

.87 

.77 

.67 

<VUt/te( 8,500)  (d-2)* 

8.74 

.09 

.87 

.87 

.64 

.48 

a>UW.e(8,500)(d-3) 

8.75 

.09 

.87 

.88 

.62 

.49 

OM.  (8,500) 

8.47 

.08 

.81 

.95 

na** 

na 

a^W.e(16,250(d-2) 

8.63 

.09 

.84 

.87 

.70 

.58 

a/UW.e( 16,250)  (d-3) 

8.74 

.10 

.86 

.98 

.75 

.65 

<we(16,250) 

8.20 

.08 

.72 

1.23 

na 

na 

/te(8000) 

8.99 

.07 

.86 

.50 

.71 

.72 

aAA/ie(8,1000)  (d-2) 

9.00 

.93 

.52 

.48 

.36 

OJiOh.l( 8,1000)  (d-3) 

9.02 

.07 

.93 

.53 

.49 

.34 

a/te(8,1000) 

8.83 

.07 

.88 

.49 

na 

na 

OAiiA£(16,500)(d«2) 

8.94 

.88 

.49 

.56 

.47 

<MUWe(16,500)(d-3) 

8.99 

.87 

.52 

.62  . 

.58 

a/te(16,500) 

8.66 

.06 

.75 

.53 

na 

na 

16000) 

9.04 

MM 

.92 

.22 

.64 

.59 

OAaAZi 8,2000) (d-2) 

9.05 

MM: 

.93 

.22 

.50 

.43 

<VUXAe.(8t2000)  (d-3) 

9.05 

.05 

.89 

.22 

.53 

.48 

aie(8,200) 

8.96 

.05 

.88 

.20 

na 

na 

0404e(16, 1000) (d-2) 

9.04 

.05 

.93 

.22 

.47 

.42 

CUUVIL(.  16, 1000)  (d-3) 

9.06 

.05 

.89 

.24 

.53 

.48 

a*e(16,i000) 

8.88 

.04 

.87 

.21 

na 

na 

<VUWe(32,500)(d-2) 

9.00 

.05 

.91 

.20 

.54 

.49 

O'iaAe(32,500)  (d-3) 

9.04 

.05 

.91 

.23 

.59 

.56 

0Ae( 32,500) 

8.71 

.04 

.79 

.26 

na 

na 

£ 

d  -  degree  of  polynomial  fitted  in  regression 


**  . , 
na:not  available 


subsection  lengths  {n^}  ■  n/10,  [n/7j,  n/5,n/4,  [n/3],  n/2,  n 
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Table  3 


Simulation  Results  for  Estimates  of  E(W)  ■  12.77  In  M/G/l  Queue 

. . . . . . . . .  „ 

with  p  -  0.9  (X  *  1.0,  y^  -  0.5,  Vj  “  2.0,  p  *  0.33).  Here  R  -  200. 


Estimate 

_ 

w 

iHl 

90% 

Coverage 

MSE 

Correlation 

(W,S(W)) 

Spearman 

Correlation 

(W.S(W)) 

*e(5000) 

.09 

.810 

1.73 

.71 

.77 

<yuwe(io,500)(d-2)* 

1 

.09 

.865 

1.78 

.61 

.60 

a/UUte(10,500)  (d-3) 

12.61 

.10 

.875 

1.89 

.62 

.60 

a/U><  10, 500) 

12.03 

.08 

.770 

2.02 

na** 

na 

ioooo) 

12.64 

.07 

.850 

mm 

.70 

.71 

a*a/te(l0,1000)  (d-2) 

12.66 

.880 

Ba 

.50 

.45 

a/uvte(io,iooo)(d«3) 

12.72 

.07 

.910 

Ba 

.52 

.46 

a/ie(l0,1000) 

12.39 

.860 

na 

na 

cuuv te(20,500)  (d-2) 

12.51 

1 

.825 

Ba 

.60 

.58 

<VtO/te(20,500) (d-3) 

12.60 

.825 

.65 

.63 

<Wfc(20,500) 

12.00 

l 

.700 

ns 

na 

na 

/ie(  20000) 

12.78 

.05 

.855 

.51 

.72 

.73 

<VUVLt( 5,4000) (d-2) 

12.80 

.05 

.885 

.51 

.39 

.36 

a/UUte(5,4000)  (d-3) 

.05 

.880 

.53 

.37 

.35 

a*e(5,4000) 

12.71 

.05 

.870 

.50 

na 

na 

<VlflAe(10,2000) (d-2) 

12.78 

.05 

.880 

.50 

.49 

.49 

ata’ie(10,2000)  (d-3) 

12.81 

.05 

.880 

.53 

.52 

.52 

<wte(io,2000) 

12.63 

.05 

.820 

.50 

na 

na 

<VUVLZ{ 20, 1000)  (d-2) 

12.73 

.05 

.860 

.52 

.52 

.49 

OAjtVUi 20, 1000)  (d-3) 

12.79 

.05 

.855 

.56 

.55 

.52 

044(20,1000) 

12.47 

.05 

.815 

.59 

na 

na 

d  «  dgree  of  polynomial  fitted  In  regression 


na:  not  available 
for  n  ■  500, 
n  -  1000, 
n  -  2000 
n  -  4000 


-  50,  71,  100,  125,  166,  250,  500 

{i^}  -  50,  100,  125,  166,  200,  250,  333,  500,  1000 

-  100,  166,  200,  250,  400,  500,  666,  1000,  2000 
{i^}  -  400,  571,  800,  1000,  1333,  2000,  4000 
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Table  4 


Summary  of  Notation 


*fc(p) 


— regenerative  estimate  based  on  p  cycles; 


— realization  of  /ifc(p)  on  the  1th  subsection  out  of  the  m^  sub¬ 
sections  In  a  replication. 

— average  regenerative  estimate  from  subsections  of  length  m^; 

m. 

-  I  -*e(J,n.)/m.. 
j-1  K 

<M.e(j,nijt,njt)  —realization  of  a^feCm^.n^)  on  the  jth  replication. 

a/ietmm^.n^)  —overall  average  regenerative  estimate  for  mm^  subsections  of 

m 

length  -  l  0/te(j,mk,nk)/m. 
j*l 


-ted,  p) 


oteCn^,^) 


— sample  variance  of  ate(mmjt>n^); 

®  — —  2  1/2 

■  {  l  (ote(j  .m^.n^)  -  OAeCmm^.nj^))  /m(m-l)}  '  . 


tate(n)  — regression  adjusted  regenerative  estimate  using  n  cycles. 

-tate(j,n)  — realization  of  4ate(n)  on  the  jth  replication  of  a  section  of 

length  n. 


otate( m,n) 


MO*e(m»n) 


— average  regression  adjusted  regenerative  estimate  from  m  sections 
m 

of  length  n;  ■  \  -tate(j,n)/m. 

J-1 

— sample  standard  deviation  of  ata*e(m,n); 

•  {  \  (*ate(j,n)  -  atat^m.nJ^/nKm-l)}1^. 


J-1 
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